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Abstract — One source of disturbance in a pulsed T-ray signal 
is attributed to ambient water vapor. Water molecules in the 
gas phase selectively absorb T-rays at discrete frequencies cor- 
responding to their molecular rotational transitions. This results 
in prominent resonances spread over the T-ray spectrum, and in 
the time domain the T-ray signal is observed as fluctuations after 
the main pulse. These effects are generally undesired, since they 
may mask critical spectroscopic data. So, ambient water vapor 
is commonly removed from the T-ray path by using a closed 
chamber during the measurement. Yet, in some applications 
a closed chamber is not applicable. This situation, therefore, 
motivates the need for another method to reduce these unwanted 
artifacts. This paper presents a study on a computational means 
to address the problem. Initially, a complex frequency response 
of water vapor is modeled from a spectroscopic catalog. Using a 
deconvolution technique, together with fine tuning of the strength 
of each resonance, parts of the water-vapor response are removed 
from a measured T-ray signal, with minimal signal distortion. 

Index Terms — T-rays, Terahertz, THz-TDS, water vapor, ro- 
tational transitions, removal 



I. Introduction 

Terahertz time-domain spectroscopy (THz-TDS) has re- 
ceived much attention from researchers due to its outstanding 
performance [1], [2]. With the assistance of ultrafast femtosec- 
ond laser technology, a T-ray system generates and detects 
a short pulse, and thus the corresponding frequency range 
spanning from a few hundred gigahertz to a few terahertz or 
more. This frequency range has a relatively high black-body 
radiation background, prohibiting a high-SNR detection with 
a conventional detector. However, this problem is suppressed 
by adopting a coherent detection scheme, lifting the SNR to 
as high as 60 dB [3]. 

In the T-ray frequency region, many polar gases of general 
interest possess unique rotational transition energies, which 
give rise to spectral resonances [4]. Because of these unique 
fingerprints, THz-TDS proves useful for gas classification and 
recognition. However, this property, on the other hand, affects 
the THz-TDS capability in an open air setting, in which water 
vapor is ubiquitous. 
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Water vapor, the third most abundant gas in nature [5], 
is known to have numerous rotational resonances in the T- 
ray band [6]. THz-TDS of a sample, in open air, therefore 
unavoidably results in a combination of the sample's spectral 
features and water vapor resonances in the frequency domain. 
In the time domain, this results in field fluctuations after the 
main T-ray pulse. Mostly, these effects are undesirable, since 
they can obscure spectroscopic results of interest. 

The water-vapor effects are removable during the mea- 
surement by purging an enclosed T-ray path with dry air 
or such a non-polar gas as nitrogen, which does not have 
transition energy levels in the T-ray regime [7] — alternatively 
a vacuum is sometimes used. In some applications, it is not 
always possible to enclose the entire T-ray beam path, for 
example, in applications where stand-off detection is required 
[8]. The effects can be partially canceled out by calculation, 
if two measurements, one for the sample and the other for 
the reference, are available. Again, not all situations can 
provide both measurements. Performing narrowband T-ray 
sensing in an atmospheric transmission window could avoid 
the effects [9], at the expense of a full spectral fingerprint. 
As a completely different approach, a numerical method to 
alleviate the water vapor effects on the data would therefore 
be beneficial. 

Knowing the resonances' characteristics in the T-ray fre- 
quency range allows a numerical estimation of the complex 
response of water vapor. Theoretically, the estimated complex 
response can be deconvolved directly from a received T- 
ray signal. But pragmatically, fitting of a modeled complex 
response to a measured response is complicated by many 
factors, e.g. the limited dynamic range [10], limited frequency 
resolution [11], and measurement uncertainties [12]. Besides, 
the exact model requires precise knowledge of geometric and 
atmospheric conditions during the measurement. Moreover, if 
direct deconvolution is carried out, there is no measure to 
validate the result. 

Fine-tuning the strengths of complex resonances based on a 
brute-force search is introduced in this work. Each resonance 
is tuned in magnitude within a predefined range and then 
deconvolved from a measured signal. A tuning criterion is 
met when the ratio of the fluctuation energy to the main pulse 
energy of an adjusted signal is minimized. Repeatedly tuning 
the strength line-by-line ultimately results in the mitigation of 
water- vapor-induced effects. Constrained by the condition of 
minimum fluctuation energy, the algorithm always provides 
the optimum results. 

The proposed algorithm based on signal processing has 
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Fig. 1. Effects of water vapor on a T-ray pulse and its spectrum, (above) T- 
ray signals recorded in nitrogen and water-vapor atmospheres. The sampling 
interval is 0.067 ps, and the total time duration is 136 ps. (below) Their 
corresponding spectra with the spectral resolution of 7.32 GHz. For this 
measurement the noise floor marks the cutoff frequency at approximately 
4.0 THz. 



merits in that (i) its generality allows application to T-ray 
signals with reduced effect on samples' temporal and spectral 
features, and (ii) the exact measurement conditions, including 
humidity, temperature, pressure, and propagation length, are 
not crucial in order to remove the water-vapor response. 

This paper is organized as follows: Section [TT] elaborates 
the effects of water vapor on T-ray signals and spectra. Sec- 
tion|IIlJ the model characteristics of each complex resonance is 
determined and validated with a measurement. In Section |IV] 
the algorithm for removal of water-vapor effects is introduced. 
This algorithm is verified by experiments in Section [V] fol- 
lowed by discussions and a conclusion in Section [VI] 

II. Effects of water vapor on pulsed T-ray signal 

Water vapor contains many modes and the details of its 
interaction with an electromagnetic wave depends on the fre- 
quency. Water vapor exhibits pure rotational modes of energy 
transitions spanning the millimeter wave to mid-infrared or 
approximately 3 GHz to 19.5 THz [13], [14]. In the infrared 
region, both pure rotational line transitions and rotation- 
vibration bands are noticeable [15]. These transitions cause 
absorption and re-emission of wave energy in a number of 
narrow frequency bands, unique to water molecules. 

Fig. [T] shows T-ray signals and spectra measured by THz- 
TDS under nitrogen and water-vapor atmospheres at room 
temperature and pressure. In the time domain, the signal 
recorded in the presence of water vapor undergoes long 
fluctuations after the main pulse. This is due to the energy 
re-emission by rotational transitions of water molecules. In 
the frequency domain, the water vapor causes several sharp 
resonances at discrete frequencies, as a result of quantised 
rotational transition energies. In terms of Fourier theory, 
the pair of time-domain fluctuations and frequency-domain 
resonances is based on the principle that a sharp feature in 
one domain is related to a broad feature in the other domain. 



In addition to the fluctuation and the resonance effects, the 
analysis shows that the T-ray energy loss by water absorption 
calculated between 0.0 and 4.0 THz is as high as 10% for the 
spectrum in Fig.[T] The loss is likely due to the non-directional 
energy re-emission of water molecules. In the time domain, the 
ratio of the main pulse energy to the tail (fluctuation) energy 
is calculated for both of the time-domain signals (for more 



details on the calculation, see Section IV-C I. The energy ratio 



for the nitrogen measurement is 429.98, whereas the ratio for 
the water vapor measurement is 18.82. This measure will be 
used later on to construct a criterion for fluctuation removal. 

III. Model of water vapor absorption and 

DISPERSION 

It is required that the water-vapor response be modeled 
appropriately before a further step in water-vapor removal is 
taken. Modeling of absorption and dispersion for a rotational 
resonance involves selecting a suitable lineshape and calculat- 
ing the line position, line strength, and linewidth. Fortunately, 
the line position and line strength are, only to a small extent, 
dependent on atmospheric conditions, i.e. the line position is 
constant at low pressures [16] and the peak absorption is not 
a function of the pressure [16], [17]. In fact, the line position 
and line strength have been well parameterized in many 
publications. Thus these two parameters are readily available, 
leaving only the linewidth and lineshape to be determined. In 
the following subsections these parameters are discussed in 
detail relative to water molecule. 

A. Line strength and line position 

Water vapor rotational resonances in the T-ray regime have 
been extensively measured and analyzed via either conven- 
tional FTIR spectroscopy [18], [19], THz-TDS [6], [20], or 
other techniques [21]-[23] — the extensive study is probably 
due to the abundance of water vapor, which involves many 
physical processes in nature. Subsequently, the computer- 
accessible spectroscopic parameters of water vapor, including 
the line position, line strengths and linewidth, are available 
from many research groups. Here in this work the catalog 
published by the JPL group [24] is adopted. As this catalog has 
been regularly updated since before 1985 [25], it is unlikely 
that the high-intensity lines of water vapor are absent from the 
list. 

According to the JPL catalog, the line strength is reported in 
terms of the integrated line intensity at 300 K, 7 Q (300 K), and 
is scalable with the temperature, T. Despite the temperature 
dependency, the line strength is independent of pressure [16]. 
Also included in the catalog is the line position. The shifting 
of the position is an order of magnitude lower than the line 
broadening [26], and thus the position is presumably pressure- 
and temperature-dependent [16]. 

B. Linewidth 

A rotational absorption resonance is possibly broadened 

by a number of factors, for instance, self or foreign-gas 

collisions, molecule-wall collisions, the Doppler effect, and 
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natural lifetime broadening [16], [27]. Among these factors, at 
a standard pressure the collisional broadening is predominant 
in determining the spectral linewidth [17], [28]. 

More specifically, the collision broadened linewidth is de- 
pendent on pressure [27] and temperature [15], [20], [26]. This 
dependency yields the Benedict and Kaplan relationship [13], 
[29], 

- - ^ i£) {ff • 

where Awrj is the FWHM of the line at pressure po and tem- 
perature T , Aoj is the FWHM at pressure p and temperature 
T, and m is the temperature index, varying between 0.5 and 
1.0. An index of m = 0.68 can be used in the calculation [30]. 

Determined from the experiment, the average FWHM, 
Awo/2ir, of water vapor resonance at ambient temperature 
and standard pressure is 6 GHz. This value is in agreement 
with an estimate of 10 MHz per Torr, or 7.6 GHz at 1 atm [27], 
[31] within the limit of measurement uncertainty. 
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Fig. 2. Model and measurement of H2O absorption. The lower graph is a 
continuation of the upper graph, as indicated by the arrowheads. The measured 
and modeled curves match well in the range of 0.0 to 1.6 THz, but match 
poorly in the higher frequency range, where the absorption strengths become 
higher than the dynamic range of the system. The maximum measurable 
absorption coefficient of the system, represented by the dashed line, is 
calculated using the method of Jepsen and Fischer [10]. 



C. Lineshape 

The rotational absorption lineshape can be modeled either 
by Lorentz [32], Vleck-Weisskopf [33], or Gross [13], [34] 
profiles. The choice depends on the frequency range of interest 
and the nature of line broadening. Notwithstanding, no signifi- 
cant difference is observable among these profiles at or nearby 
a T-ray resonance, with the linewidth nearly three orders of 
magnitude lower than the resonance frequency. In this paper, 
a Lorentzian profile is adopted in modeling the water vapor 
response in the T-ray frequency range. 

A Lorentz absorption and dispersion profiles are given 
by [13], respectively, 
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where uj a denotes the a th transition frequency, and Aui a 
denotes the HWHM of the profile at that frequency. The 
absorption and dispersion profiles are linked together via 
Kramers-Kronig relation. 

D. Ensemble of rotational transition resonances 

This subsection gives a complete model, in terms of optical 
constants, as a result of several rotational transition resonances 
in the frequency range of interest. The model is a function of 
the lineshape, linewidth, line strength, and line position, which 
are discussed in the previous subsections. Later on, the model 
is compared to the measurement to verify the accuracy. 

The model extinction coefficient and index of refraction are 
summation of an offset and a set of Lorentzian profiles, or, 
respectively, 

k(u) = k(0) + ]T M , (3a) 

a 

n(w)-l = n(0)l^m„KH-l], (3b) 



where the line strength m a is proportional to the integrated 
line intensity, I a , given in the JPL catalog, divided by the 
line position, u> a , or, m a cx l a juj a . The offset is partially a 
result of electronic and molecular vibrational resonances at 
high frequencies [13], and also a contribution from the tails 
of other rotational resonances unaccounted for by the model. 

Typically, at the low frequency range, from DC to a few 
hundred gigahertz, THz-TDS cannot produce sufficient energy 
to overcome the noise. As a result, the resolved optical 
constants in this frequency range are unreliable. A phase 
extrapolation technique is usually introduced to correct the 
unwrapping process, forcing the phase to start at zero [35]. 
Thus, it is not necessary in the model to consider the index 
offset, n(0), which is derived from the phase. 

The absorption of water vapor, shown in Fig. [2] is deter- 
mined from the measured data (Fig. [TJ and the Lorentzian 
model in the frequency range between 0.0 and 4.0 THz. It 
can be clearly seen that below 1.6 THz the model closely 
resembles the measurement. This match is possible since the 
resonances in this low frequency range have their strengths 
beneath the maximum absorption coefficient measurable by 
the system. However, as the frequency goes beyond 1.6 THz, 
the model cannot track the measurement due to the dynamic 
range of the system. In this situation, the system can no longer 
measure the absorption coefficient correctly, and the clipped 
absorption peaks are obvious, in particular, in the range from 
2.0 to 4.0 THz. Furthermore, the clipped lines become less 
distinct, when the T-ray power reaches the noise floor beyond 
3.0 THz. 



E. Continuum absorption 

The continuum absorption, unaccounted for so far, is defined 
as the excess measured absorption unable to be quantified 
by the resonance spectrum. The mechanism underpinning the 
continuum has not been fully understood [36], and most of the 
mathematical models are fitted to experimental observations. 
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It is known that the continuum absorption is pressure and 
frequency dependent. A continuum for H2O collisions, derived 
from an experiment in the T-ray frequency range, between 0.4 
and 1.83 THz, is estimated at 4.22 x 1(T 8 (dB/km)/(hPa GHz) 2 
[36]. A normal atmospheric measurement in a relatively short 
path length is insignificantly influenced by the continuum 
absorption [36], [37]. 

IV. Removal of H 2 response by strength tuning 

This section proposes a numerical means for removal of 
water-vapor effects, as an alternative to the conventional 
'sealed chamber' method. Subsection II V-AI establishes the T- 
ray spectroscopic system, and points out the difficulties in us- 
ing a direct deconvolution to remove the water vapor response. 



The strength-tuning algorithm is introduced in Subsection IV 
[B] followed by the fluctuation ratio as a tuning criteria in 
Subsection IV-C Discussion on the algorithm's generality can 
be found in Subsection IIV-DI 



A. Water vapor as a black box 

During a THz-TDS measurement, an emitted T-ray signal 
evolves based on several factors, e.g. the sample response, 
the water vapor response, the instrument response, and the 
noise. These can be modelled as a system, as illustrated in 
Fig. [3] Given that x(t) denotes the input T-ray pulse, which 
is immediately deployed from a transmitter, s(t) the impulse 
response of a sample, w(t) the impulse response of water 
vapor, u(t) the impulse response of the instrument, n(t) the 
noise, and r(t) the time windowing, the received pulse, y(t), 
can be described as 

y(t) = [ x (t) * s(t) * w(t) * u(t)] ■ r(t) 

+[n(t) *«(*)]• r(t) , (4) 

where * is the convolution operator. The noise is simply 
additive and can thus be treated as an extra term in Equation [4] 
Via Fourier transform, the above equation in the frequency 
domain is 

Y(lo) = [X(u) ■ S(u) ■ W{uj) ■ U(lu)] * R(lo) + N(u) , (5) 

where Y(w), X{w), S{uS), W(w), U(w), R(u), and N(w) 
represent complex frequency responses, as the Fourier pairs 
of y(t), x(t), s(t), w(t), u(t), r(t), and [n(t) * u(t)] ■ r(t) 
respectively. 
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Fig. 3. Model of a spectroscopic system. The system response is the 
combination of the frequency responses of the sample, water vapor, and 
instrument, denoted by S(ui), W(ui), and U(u>), respectively. The symbols, 
X(u>), Y[ui), N(lo), and R(u>) represent input T-ray spectrum, measured 
spectrum, noise, and window (or limited recording duration) responses, 
respectively. The water-vapor response in the dotted box needs to be removed. 



The water-vapor frequency response can be expressed by 

p(T) luL 



W(lo) 



exp 



-jn y 



Vo(T) 



(6) 



Here, L denotes the T-ray propagation length in free space, 
excluding the space occupied by sample(s). The water vapor's 
complex index of refraction, rj vap = n vap — jK va p, contains 
the refractive index n vap and the extinction coefficient K vap , 
which are modeled in Equation [3] The water- vapor response 
is also determined by the humidity [38], where p(T) is the 
vapor density at temperature T, po(T) is the saturation vapor 
density at the same temperature. 

In order to remove the effects of water vapor from a 
measurement, the water-vapor response, W(cj), in Equation [5] 
must be replaced by the vacuum response, V{iS), with the 
same propagation length, 



exp [-jwL/ c] 



(7) 



From Equation [5] the noisy response can be modelled in two 
parts, 

f [X(u)-S(u)-W(u))-U(u)]*R(u}) 
Y(uj)^1 ^ ,if \X-S-W\>\N\ (8) 

[ N(u>) , otherwise. 

Assuming the noise level is unity and the model water-vapor 
response matches well with the measurement, performing the 
direct deconvolution would yield 

( [X(lo) ■ S(uj) ■ V(u) ■ U(lu)] * R(w) 
Y(lo)^1 ,if \X-S-W\>1 (9) 

[ V{ijj) /W(lj) , otherwise. 

Therefore, removal of the water vapor response fails at 
\X-S-W\< \N\ or D(w) < \S ■ W\~\ where D(w) = 
|X/iV| is the system's dynamic range. This usually occurs 
at the peaks of absorption resonances, as shown in Fig. [2] 

Consider the water vapor response closely. It can be sepa- 
rated into several components, and the separation is formulated 
as 



W{u) 



V(w)W (u)Wi(u) . . .W a (w) 



(10) 



Each component corresponds to a complex resonance, and has 
its frequency response derived from Equations [3] and [6] or 



W a (w) 



exp [-j m a (n a - jK a )aj/c] 



(11) 



The line strength, m a , now incorporates the humidity ratio, 
p(T) / po(T), and the propagation length, L. It can be inferred 



from Equations 10 and 11 that all components or resonances 
can be removed from the measured spectrum separately and 
independently. Hereby, the deconvolution at the noisy parts of 
the spectrum can be avoided, and also it is possible to tune 
m a of each resonance to the measurement without accurate 
information of the atmospheric and geometric parameters. 
In order to achieve this parameter relaxation, somehow, a 
criterion indicating the fit of tuning is essential. 
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B. Fluctuation- removal algorithm 

The fluctuation-removal algorithm is shown in Fig. |4] The 
counters a and b for the line position and line strength are 
reset. The H2O spectral line parameters, including intrinsic 
line positions and strengths, in the frequency range of inter- 
est, are then fetched from an existing database. A complex 
resonance profile, n(ui) — jn(uj), of a rotational transition 
at uj a =o is modeled according to the Lorentzian profile in 
Equation [2] The modeled profile multiplied by the initial 
strength factor rn a= o,6=o is then temporarily deconvolved 
from the measured complex response, Y(u>). The deconvolved 
time-domain signal, y(t), is estimated for its fluctuation ratio 
(see Subsection IV-C[ ). The procedure repeats to find the 
fluctuation ratio at other predefined line strengths m a= oy, 
6e{l,2,3,...}. 

Once the tuning range, m a= o.b, is covered, the algorithm 
picks up an optimal strength within m a= o^, which gives the 
minimum fluctuation ratio, if any, and permanently removes 
that optimal complex resonance from the measured signal. 
The tuning procedure starts again, but with the next transition 
resonance, ui a =i, and is repeated until all resonances are opti- 
mized and possibly removed. It should be noted that a strong 
resonance, which suffers from noise i.e. \X ■ S ■ W\ < \N\ 
at and around the peak of resonance, will be avoided by the 
tuning process, because performing deconvolution of such an 
ill-defined resonance will give rise to a high fluctuation ratio. 
Since there might be dependencies among the resonances, with 
regard to time-domain fluctuations, the whole process starts 
again until the fluctuation ratio no longer decreases. 

It is advisable not to perform zero padding prior to Fourier 
transform in the removal process, as the interpolated spectrum 
does not exactly reflect the reality — this is especially the case 
for points at resonances that have high variability. Otherwise, 
discrepancy between the model and the spectrum would cause 
a large remnant in the spectrum after deconvolution. Also note 
that to speed up the algorithm, the resonances that have a 
strength lower than the amplitude uncertainty can be skipped. 
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Fig. 4. Fluctuation-removal algorithm. Each resonance in the T-ray frequency 
range is modeled and fine-tuned in its strength with the criterion of minimum 
fluctuation, and then the optimal resonance is removed from the measured 
complex response. The procedure is repeated until all resonances within the 
frequency range of interest are investigated, and until the fluctuation ratio no 
longer decreases. 



C. Fluctuation ratio 

A criterion related to the quality of a T-ray signal is 
necessary in selecting the optimal resonance line strengths. 
The 'quality' here is defined as having a high pulse energy 
with low fluctuations in the time domain, which implies the 
absence of water-vapor resonances in the frequency domain. 

In order to quantify the quality one must be able to evaluate 
the total energies of the main pulse and of the fluctuations. 
One potential way is to window the time-domain signal with 
a Gaussian profile to obtain the amplitude at a desired portion. 
A Gaussian window is selected, because a T-ray main pulse 
is essentially derived from an optical pumping pulse, which 
takes on a Gaussian pulse profile [39]. 

Fig. [5ja) shows a Gaussian window, g a (t), overlapping the 
T-ray signal, y(t). A Gaussian window with appropriate width 
and position would eliminate the fluctuating tail of a signal, as 
shown in Fig. |5Jb). On the other hand, the complement of a 
Gaussian window could also be used to remove the main pulse, 
as in Fig. [5jc). With the assistance of a Gaussian window, 



the fluctuation energy normalised by the main pulse's total 
energy — coined the fluctuation ratio — can be formulated as 



F(y,g ;o) 



J t [y(t) ■ 9At)] 2 dt 



(12) 



The integration is carried out over the time duration of a 
recorded T-ray signal. A Gaussian window is given by 

g a {t) = exp[-(i - t a ) 2 /2a 2 } , (13) 

where to is the peak position of the T-ray signal, y(t), and a 
multiplied by 2V 2 In 2 is the FWHM of a Gaussian window. 

D. Generality of the algorithm 

The algorithm is plausibly general — in the sense that it can 
be applied to any T-ray signal with minimal disturbance of 
desired signal features. This generality is due to the following 
facts: 

First of all, a T-ray spectrum is altered by the algorithm 
only the frequencies at which water-vapor absorption lines 
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(c) T-ray signal weighted by the Gaussian window complement 

Fig. 5. Gaussian window applied to the T-ray signal. A T-ray signal weighted 
by a Gaussian window, which has its peak position set exactly at the main 
pulse peak, yields the main pulse with suppressed fluctuations. A T-ray pulse 
weighted by the Gaussian window complement yields only fluctuations. The 
FWHM of the Gaussian window used here is 8.33 ps. 



are situated. The spectral resonances of other polar gases, 
occupying narrow frequency bands and usually not overlap- 
ping with water resonances, are not disturbed. In addition, the 
broad spectral features of a time-domain transient is unlikely 
to be significantly affected by narrow water line removal. The 
reflections, for example, which exhibit fringes over a broad 
frequency range, are not affected by narrow band resonance 
removal. 

Second, unrelated fluctuations are ignored. Because the 
algorithm senses the fluctuation change, caused by strength 
tuning, rather than the fluctuation itself, any unrelated fluctu- 
ations cannot deceive the algorithm. 

Last, the causality of a signal is preserved. The absorption 
and dispersion profiles, used in the model of water vapor 
response, exactly comply with the Kramers -Kronig relation. 
Deconvolution of the measured spectrum by this causal re- 
sponse would yield the causal spectrum. 

V. Results 

A. Free-path measurement 

The first T-ray signal, subject to the water-vapor removal 
algorithm, is measured in free space without the presence of 
any material. The signal has the temporal resolution of 66.7 fs 
and the total duration of 136.63 ps, providing the spectrum 
with the spectral resolution of 7.3 GHz. 

The strength-tuning algorithm is carried out with a set 
of complex resonances in the frequency range between 
and 4 THz, where the T-ray magnitude is relatively high. 
The sequence, in which the resonances are interrogated and 
removed, follows the numerical order, i.e. from low to high 
frequencies. Only the resonances that have the strength higher 
than 0.01 of the maximum strength are inspected. The FWHM 
of the Gaussian window is 3 ps, and the iteration of the 
algorithm is set to five. 
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Fig. 6. T-ray signals measured in free path. The sampling interval is 66.7 fs 
with the total duration of 136.63 ps (only the first 80 ps is shown here). The 
inset shows a zoom-in of signals between 5 and 10 ps. The fluctuation ratios 
of the N2 -atmosphere, F^O-atmosphere, and F^O-removal measurements 
are 2.2, 9.2, and 3.3, respectively. The fluctuation after the main pulse is 
remarkably reduced. 



As shown in Fig. [6| the algorithm can significantly reduce 
the time-domain fluctuations, which are previously located 
immediately after the main pulse, i.e. after 10 ps. The main 
pulses in the inset clearly illustrate the similarity between 
the fluctuation-free signal and the processed signal, despite 
a temporal shift which is not accounted for by the algorithm. 
In Fig.|7J the spectrum result demonstrates success of the algo- 
rithm in removing, or reducing the strength of, the resonances 
in the frequency range from to 2.5 THz. However, some 
resonances still persist, where they are spectrally congested or 
severely disturbed by noise. 

B. DL- Phenylalanine measurement 

The second experiment is performed with a signal measured 
from a DL-Phenylalanine sample [40] at 220 K. The cryo- 
stat containing the sample is evacuated, so no water vapor 



nitrogen atmosphere 

— water-vapor atmosphere 
with fluctuation removal 
water-vapor atmosphere 




1.5 2 2.5 
frequency (THz) 

Fig. 7. T-ray spectra measured in free path. The spectral resolution is 
7.3 GHz. Up to 2.5 THz most of the water resonances are removed, or at 
least reduced in the strength, by the proposed algorithm. The magnitudes are 
offset for clarity. 
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Fig. 8. T-ray signals measured with DL-Phenylalanine sample in place. 
The sampling interval is 66.7 fs with the total duration of 68.38 ps . The 
reflections at 32 and 63 ps due to the cryostat's windows can be recovered 
from fluctuation by the algorithm. The fluctuation ratios of the N2 -atmosphere, 
H20-atmosphere, and H20-removal measurements are 2.6, 7.9, and 3.8, 
respectively. 



at cryogenic temperature involves in the measurement. The 
interacting water vapor is from normal air, surrounding the 
cryostat. The total scan is 68.3 ps at the step size of 66.7 fs. 
The spectral resolution is 14.6 GHz. 

The recorded signal particularly contains the reflections, 
caused by the cryostat's windows, made from a cyclo-olefin 
copolymer (Topas). Use of such a signal demonstrates how 
well the proposed algorithm can deal with water-vapor fluctu- 
ations without disturbing unrelated reflection features. 

The parameters set for the algorithm are similar to the 
previous test, except for the interrogated frequency range 
which is in between and 2 THz. 

The processed signal is shown in Fig. [SJ along with the 
original. Again, it is obvious that the fluctuations located after 
the main pulse, beyond 15 ps, are remarkably reduced, in spite 
of an introduction of small oscillation at the beginning. The 
reflection at 63 ps (point B) is not disturbed by the algorithm, 
and, surprisingly, the reflection at 32 ps (point A), which was 
initially buried in fluctuations, is now recovered. In Fig. [9] the 
algorithm can remove most of the H2O resonances. Still a few 
resonances cannot be completely removed. These persistent 
resonances are around 1.1 and 1.7 THz, exactly the same 
positions as the unremoved resonances in the previous case 
(see Fig. [7]). The recovery of reflections suggests a promising 
application of the algorithm in T-ray range finding. 

C. Evaluation of the algorithm 's performance 

Table [I] shows the fluctuation ratio and mean squared error 
(MSE) of the signals and the total energy of the spectra, in- 
vestigated in the previous subsections. Two additional results, 
not shown in this paper, are tested with lactose's reference and 
sample signals. The fluctuation ratio helps measure a change 
in the fluctuation energy and main pulse energy, before and 
after the signals are processed. The MSE evaluates the fitness 
of the improved signals to the N2 -atmosphere measurement. 



„ nitrogen atmosphere 

— water-vapor atmosphere 
with fluctuation removal 

— water-vapor atmosphere 




1.5 2 2.5 3 
frequency (THz) 

Fig. 9. T-ray spectra measured with DL-Phenylalanine sample in place. The 
spectral resolution is 14.6 GHz. From to 2 THz most of the resonances are 
removed, or at least reduced in the strength, by the proposed algorithm. The 
magnitudes are offset for clarity. 

TABLE I 

Quantitative evaluation of the algorithm's performance. The 
notation h2o, and h2o stand for measurements in water vapor 
and water vapor with fluctuation removal, respectively. 



Signal 


Fluctuation ratio 


MSE (%) 


Energy (a.u.) 


H 2 


H 2 


H 2 


h 2 o 


h 2 o h 2 o 


Free path 


4.2 


1.5 


0.076 


0.047 


0.90 0.96 


Phenylalanine 


3.0 


1.5 


0.13 


0.085 


0.90 0.95 


Lactose reference 


16.1 


2.0 


0.31 


0.19 


0.90 0.98 


Lactose sample 


4.8 


1.7 


0.31 


0.21 


0.93 0.99 



The MSE is given by 

MSE = 



7T) * * 



Vin) 2 , 



(14) 



where y is a signal measured in the nitrogen atmosphere and 
y is a compared signal. The error is summed and averaged 
over m temporal points. The total spectrum energy shows the 
recovery of the energy from resonance absorptions. Note that 
the fluctuation ratio and the spectral energy are normalised to 
those of the signal measured in nitrogen atmosphere. 

From the table the fluctuation ratio indicates that the fluctu- 
ation energy is greatly reduced in all cases, once the removal 
algorithm is implemented with the signal measured in the H2O 
path. However, the MSE still reflects a considerable difference 
between the processed signals and the signals measured in 
a nitrogen-filled chamber, in spite of the great improvement 
visualized in the previous subsections. This would probably 
be due to the whole signal shift in the time domain, caused 
by a constant refractive index difference unaccounted for by 
the algorithm. 

The total spectral energy in Table [I] reveals the algorithm is 
successful in recovering a part of the signal's energy, which 
is absorbed by water vapor resonances during propagation. 
However, not all the energy is recoverable. With regard to the 
processed spectra shown in the previous subsections, the algo- 
rithm cannot completely remove the water vapor resonances 
in the region where the SNR is low and where the spectral 
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lines are crowded. The causes of these two cases are analyzed 
separately as follows. 

At a low SNR region, distortions in the absorption and 
dispersion resonances by noise causes deviation of the mea- 
surement from the model. This situation is demonstrated in 
Fig. [2] Applying the removal algorithm to a low SNR region 
would be effective, provided that the distortion is not severe. 
Otherwise, the large deviation inhibits the removal via a sub- 
optimal value of the fluctuation ratio. A possible means to 
alleviate the noise issue is to implement a signal enhancement 
methodology such as wavelet denoising [41] prior to the 
removal process; this might better recover water resonances 
and subsequently improve the removal efficiency. 

Crowded resonance features result in overlapping and merg- 
ing between two or more resonances. If the spectral resolution 
is too low, these blended resonances are rendered incorrectly, 
resulting in distorted shape. The algorithm would see them 
as a single stronger and wider resonance, and as a result, the 
removal does not perform well. Improving the spectral resolu- 
tion would allow better removal performance. For example, the 
free space measurement in Fig.|7]shows a better result than the 
DL-Phenylalanine measurement in Fig. [9] in the regions around 
1.1 and 1.7 THz; the former and latter measurements have the 
spectral resolutions of 7.3 GHz and 14.6 GHz, respectively. In 
addition to the spectral resolution problem, blended resonances 
cannot be represented exactly by a simple sum between two 
Lorentzian lines [42]. 

VI. Conclusion 

A numerical algorithm for removal of water- vapor effects in 
T-ray measurements is investigated in this paper. Given only a 
sample signal with no reference, theoretically, the removal can 
possibly be carried out by simple deconvolution of the model 
complex water-vapor response from the signal. However, many 
factors limit the usability of the simple deconvolution. An 
exact model of water vapor resonance requires many phys- 
ical data, including the pressure, temperature, humidity, and 
propagation length. In addition, if the noise level is sufficiently 
high such that some strong resonances are distorted, the whole 
deconvolution process cannot be performed. 

The proposed algorithm fine-tunes the strength of each 
model spectral resonance. A criterion for strength tuning is met 
when the fluctuation ratio of the processed signal reaches the 
minimum. Once an optimal resonance is attained, it is removed 
from the signal. The algorithm then proceeds to removal of the 
next resonance. This tuning scheme relaxes the requirement for 
precise information of the delicate physical conditions during 
the measurement. Furthermore, the fluctuation ratio criterion 
inhibits any fault deconvolution that might occur when the 
noise disturbs the quality of a measured signal. 

In the experiments, the algorithm produces promising re- 
sults. It can reduce a significant amount of the signal fluctua- 
tions and spectral resonances with small disturbances to other 
non-related features, such as reflections or sample-induced res- 
onances. Moreover, some features, which are initially masked 
by fluctuations, can be recovered by the algorithm. Optical 
constants extracted from improved signals have an acceptable 
quality. 



The proposed algorithm is general in the sense that any 
other polar gas response could be targeted, in principle, and 
hence removed from a measured spectrum if desired. It is not 
necessary that the parameters of a targeted gas be present in a 
spectroscopic catalog, as long as its pure response in the fre- 
quency range of interest is fully estimable. This scheme might 
have benefits in some particular situations, for instance, where 
molecular contaminant(s) are unavoidable in measurement. 

Finally, it should be noted that the algorithm is not meant to 
serve as a full replacement to the 'sealed chamber' procedure 
executed during the measurement stage. In fact, it assists the 
measurement in cases where a sealed chamber is not feasible. 
Still the algorithm cannot remove some traces of resonances 
due to the ubiquitous noise and resonance crowding. Further 
improvements of the algorithm remains an attractive challenge. 
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